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INTRODUCTION 


The  flow  in  transonic  turbomachinery  is  one  in  which  most  of  the 
difficult  complications  in  fluid  mechanics  occur  simultaneously.  The  flow  is 
in  general  transonic,  three-dimensional,  and  unsteady;  the  effects  of  viscosity 
are  felt  in  shock-wave  and  boundary- layer  phenomena  which  are,  themselves, 
three-dimensional  and  unsteady.  And  the  three-dimensionality  may  appear  in 
the  form  of  velocity  distortions  and  entropy  gradients  in  the  inlet  flow  itself. 

Despite  these  complications,  axial-flow  machines  run  with  very  attrac- 
tive levels  of  performance.  However,  increasing  demands  for  improved  performance, 
with  decreased  weight,  noise,  and  fuel  consumption , J ^ ^ continue  to  place  great 
emphasis  on  the  importance  of  understanding  these  complex  flow  problems. 


In  seeking  approaches  that  may  assist  this  understanding,  it  is 
natural  to  turn  to  methods  which  have  already  proven  their  value  in  studies  of 
the  external  flow  over  wings  and  bodies.  These  methods  can,  in  some  cases,  be 
adapted  to  the  special  circumstances  of  turbomachinery  flows;  when  this  can  be 
done,  a wide  variety  of  analytical  tools  and  a great  wealth  of  understanding 
is  brought  to  bear  on  the  problem. 

The  present  study  is  one  such  attempt.  During  the  past  decade,  the 

field  of  external  transonic  flows  has  undergone  a revolution,  centered  on  the 

use  of  numerical  methods  of  flowfield  prediction.  One  of  the  key  contributions 

that  set  these  advances  in  motion  was  the  development,  by  Krupp,  Murman,  and 
4-7 

Cole  of  relaxation  methods  that  were  capable  of  handling  the  mixed  elliptic/ 
hyperbolic  character  of  transonic  flow.  These  early  works,  all  done  in  the 
nonlinear  small-disturbance  approximation,  clarified  many  features  of  transonic 
flow,  such  as  shock-wave  and  sonic-line  locations,  and  the  relations  between 
surface  geometry  and  surface  loading.  They  also  served  as  the  necessary  external- 
flow  adjunct  for  boundary -layer  studies,  and  fo’*  predictions  of  shock-wave/ 
boundary- layer  interaction  phenomena.  The  experience  gained  in  the 
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small-disturbance  approximation  soon  led  to  advances  into  the  more  fully  non- 
linear aspects  of  the  problem,  and  into  such  further  complications  as  three- 
dimensionality. 

The  basic  purpose  of  the  present  work  is  to  apply  that  same  first 
step  to  turbomachinery  flows,  i.e.,  to  develop  a three-dimensional  relaxation 
method  for  transonic  small -disturbance  flow  through  a turbomachine  blade  row. 

The  simplicity  offered  by  the  small-disturbance  approximation  facilitates  the 
adaptation  of  external -flow  computational  methods,  and  indicates  fruitful 
directions  for  higher-order  extensions.  At  the  same  time,  the  restrictions 
implicit  in  the  small-disturbance  formulation  limit  the  range  of  validity  of 

I 

its  results  to  low  turning  angles  and  total  pressure  ratios  near  one.  The 
precise  limits  of  this  range  will  have  to  be  determined,  by  comparison  with 
experiment  and  with  the  predictions  of  more  complete  models.  In  the  meantime, 
results  found  in  the  small-disturbance  approximation  can  be  expected  to  produce 
useful  information  about  three-dimensional  shock  patterns,  boundary- layer 
behavior,  and  the  interaction  between  adjacent  supersonic  and  subsonic  zones. 

This  paper  begins  with  a summary  and  description  of  how  the  small - 
disturbance  formulation,  given  in  Part  I (Ref.  8 J , has  been  cast  into  finite- 
difference  form.  This  section  also  contains  comments  on  some  of  the  essential 
programming  steps  that  went  into  the  computer  code  developed. 

Section  II  contains  results  from  two  demonstration  calculations  that 
were  made  for  the  purpose  of  illustrating  the  computer-program  capability. 

These  results  include  examples  of  supercritical  flow  in  the  blade  passages, 
and  serve  as  a basis  for  estimating  the  computing  time  required  to  do  a given 
calculation . 

The  concluding  remarks  indicate  two  avenues  of  research  which  are 
opened  up  by  the  successful  demonstration  of  this  technique.  One  of  them  lies 
in  the  general  direction  of  further  extensions  of  the  method;  this  area  in- 
cludes such  topics  as  the  incorporation  of  the  full  nonlinearity,  and  an 
accounting  for  some  of  the  shock-wave  and  boundary- layer  phenomena  that  play 
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such  an  important  role  in  turbomachinery  performance.  The  second  direction 
consists  of  applications  of  the  present  capability.  The  availability  of  this 
method  makes  it  possible  to  conduct,  at  very  modest  cost,  parametric  design 
studies  of  the  blade  shapes  required  to  produce  a given  loading,  and  of  the 
off-design  performance  of  these  blade  rows,  all  in  the  presence  of  three- 
dimensionality  and  the  occurrence  of  mixed  subsonic/supersonic  flow.  While 
these  calculations  are  limited  to  the  range  of  validity  of  the  small-disturbance 
equations,  they  nevertheless  offer  an  immediate  method  for  learning  in  detail 
about  the  intricacies  of  this  complicated  flow  field.  There  is  much  to  be  learned 
from  such  an  application,  and  the  opportunity  need  not  be  deferred  until  after  the 
more  advanced  approximations  come  into  existence. 

i 

i 

I 

I 

I 
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I.  FINITE-DIFFERENCE  FORMULATION 


The  coordinate  system  and  rotor  geometry  used  are  shown  in  Figure  1, 
where  the  dimensionless  variables  are  defined  as 
uj  x 

2 = — , p = — , % = e - z ■,  <p  = — - (i) 

a*> 

The  velocity  components  seen  by  a blade-fixed  observer  are 

VJ x - + u , VJr  = v , = uJr  + u/  (2) 

These  dimensional  perturbation  velocities  are  related  to  the  velocity 
potential  by 

u.  3<t>  \ d<t>\  dtp  \ 


LL 

9±\  _ 

dtp  \ 

V 

dtp  \ 

dp  \ 

dp'it, 

oJ 

± d<p  ) 

t dtp  ' 

^oO 

P deLP  ' 

P d £,  / 

The  velocity  potential  satisfies  the  equation 

{’  - Mco  +P'L)  - M+  ')  M*  <Pi  1 <PiZ  t y024>„  - 2 ( I + p1)  <Pgt 
1 ] (4) 

* -ff2^  * <’V>  *,)  - o 

Subscripts  denote  derivatives  in  the  2 , p , % coordinate  system;  thus  the 

symbol  4>  stands  for  -4-^-")  . There  are  8 blades  in  the  row,  and  they  are 

7 z 3 H J/o,A 

taken  to  lie  in  the  helical  surfaces  defined  by  and  uJr  . The  axial 

projection  of  their  chord  is  a constant,  . Thus,  they  are  located  at 


062*—-  ; ± p 6 Pj  ; % 

'-*'00 


1 = 0,  1 , Z,  ...  , 8-7 
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For  purposes  of  numerical  calculation,  attention  is  restricted  to  the  region 
between  two  blades  ( 0 ^ b $ -g- ) , and  the  solution  is  made  periodic  in 
successive  intervals  of  zrc/ S>  . In  the  small-disturbance  approximation,  the 
blade-surface  boundary  conditions  are  enforced  in  the  surfaces  £ = o and 
£ ■=  2TX./B  . 

The  region  in  which  the  solution  is  to  be  found  is  now  divided  into 
a grid,  as  shown  in  Figure  2,  where  the  indices  L , K , and  N are  used  to 
number  the  points  in  the  % , 2 , and  p directions,  respectively.  Values 

of  the  potential  at  each  of  these  grid  points  are  denoted  by 

<f>  (2  , 0 , K ) = <t>(L,K,  N)  = V 

The  use  of  this  grid  makes  it  very  easy  to  enforce  the  periodicity  conditions; 

for  example,  when  solving  for  the  flow  in  the  blade-to-blade  plane  at  constant 

p , the  periodicity  condition  upstream  of  the  blade  row  (see  Equation  49  of 

Reference  8)  requires  that  <f>  (±. , K , N ) = <P  ( i-MX  , K , N)  , where  L = 1 and 

L=  Lf* IX  are  the  indices  corresponding  to  % - 0 and  % = 2tz/3  . However, 

the  skew  coordinate  system  leads  to  a misalignment  between  the  fore  Mach  cone 

through  a given  point  and  the  grid  points  in  the  ^-direction.  It  is  known 

9-11 

from  external-flow  studies  that  a misalignment  of  this  sort  will  lead  to 
an  unstable  iteration  process,  even  though  derivatives  in  the  2 -direction  at 
constant  £ are  taken  in  the  upwind  direction. 

The  remedy  for  this  problem  is  to  use  the  rotated-difference  scheme 

9 

introduced  by  Jameson.  His  method  was  developed  for  the  case  where  the  flow 
velocity  displays  large  deflections  from  one  of  the  coordinate  directions  (see 
Figure  3,  which  shows  a rectangular  coordinate  system  for  convenience). 

Jameson's  basic  contribution  was  to  replace  the  partial  differential  equation 
in  x and  by  one  in  the  streamline,  normal  coordinates  S , n . The  potential 
equation  can  then  be  expressed  locally  in  the  form 


where  the  three  dots  stand  for  lower-order  derivatives,  and  where  the  second 
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derivatives  are  related  to  their  local  counterparts  in  the  rectangular 
system  by 

tss  = J,  (?*  *«  + ***  Mj 

^ r”n  “ ^ H-j  *i  + ) 

where  q.  is  the  magnitude  of  the  local  velocity  vector  and  <jx  and  are 

its  x and  y components.  After  expressing  the  potential  equation  this  way, 
contributions  to  the  <pss  term  are  tnen  represented  by  upwind  differences, 
while  terms  contributing  to  are  evaluated  by  central  differences. 


For  the  present  case,  the  misalignment  between  the  fore  Mach  cone  and 
the  2 = constant  grid  points  is  not  caused  by  large  flow  deflections,  but  by 
the  skew  nature  of  the  grid.  The  same  general  principles  can  be  applied, 
however;  the  contributions  to  dfss  and  djnr,  can  be  identified  by  writing 
(in  dimensional  variables): 


^ <,  (lVr)Zcf,ty,  - 


W A 

+ ~rr 


where  WQX  = * (cor . Note  that  the  small -disturbance  approximation 

has  already  been  used  here,  by  neglecting  the  contribution  of  the  perturbation 

velocities  to  these  formulas.  If  these  are  now  re-written  in  dimensionless 

form,  and  transformed  to  the  2 , /£>  , t?  coordinates,  it  is  found  that 

z 

u. ) us 

= 1 4-di  > 

Thus,  it  is  clear  that  the  potential  equation,  as  displayed  in  Equation  (4), 
is  already  in  a form  suitable  for  the  application  of  Jameson's  technique:  the 

first  term  in  Equation  (4)  is  to  be  evaluated  by  upwind  or  central  differences, 
depending  on  whether  the  local  Mach  number  is  supersonic  or  subsonic,  while  all 
the  remaining  terms  are  to  be  evaluated  by  central  difference  formulas. 


+■  ( tor  ) ‘ 
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In  addition  to  these  considerations,  it  is  essential  that  the 
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difference  scheme  be  fully  conservative,  ’ in  order  to  avoid  the  spurious 
introduction  of  additional  mass  flow.  Thus,  the  potential  equation  is  re- 
written in  the  divergence  form: 

v £".*«*]■-&■{<•<"/>  hK-  ZTf  ♦,]} 

* JL  fJL  4,  . J-  4,  ll , llfl  1.  (fit,,  ■ 0 (7) 

t P M * "P  P >P 

The  grid  used  throughout  the  present  work  has  uniformly  spaced  points  in  the 
radial  and  blade-to-blade  directions , and  a nonuniform  spacing  in  the  t -direction. 
The  H -coordinate  is  taken  to  be  related  to  a transformed  coordinate  t , 
such  that 

i , . t _ dt 

di  r n ’ di  w 

A constant  grid  spacing  At  is  then  used,  giving  rise  to  a nonuniform  spacing 
in  2 . The  computer  program  as  written  contains  a specific  transformation, 

described  below,  in  one  of  the  subroutines.  Other  transformations  can  be  used 
by  making  appropriate  changes  in  this  subroutine,  whose  purpose  is  to  return 
the  function  f . 

Using  the  variable  Z , the  potential  equation  becomes: 

i f „ r o,-1  e ^ 1 

- j t\t-H.O-P  f j^\ 

n ip*'- 77^77  (f  17  If  7 Ti  17  *'  ../>■  nf 

+ 7 f ip  ' * (9) 

The  first  term  in  this  equation  is  to  be  evaluated  by  type-dependent  differences. 
In  order  to  retain  mass  conservation,  it  is  necessary  to  use  the  shock-point 
operator  at  points  where  the  flow  changes  from  supersonic  to  subsonic  values.1" 
Jameson14  has  pointed  out  a simple  way  of  doing  so,  using  a switching  function 


i-  - fL  ill 

1 p ( t.p*  it j 
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pi  which  is  defined  to  be  zero  in  subsonic  flow  and  1 in  supersonic  flow. 
Using  these  concepts,  plus  the  formulas  suggested  by  South  and  Jameson,11  the 
following  finite-difference  form  of  the  potential  equation  results: 

+ k-  <fK.*  - ~k.,)  - f.-aCk.,  - “Cj)} 


2 A £ 

- 1 

O + p1)1 

/ Ar 

P2f< 

( 

o +/01; 

/ At 

f 

rK 

l A p 

_ A/J.L 

_L\ 

L we  ( 

■f 

^.var< 

Ni  i- 

" <t>K 

4- 

■f 

- " <t>L  *'  + 

- T 

^ K- 

+ f 

WiL*l 


where  /30  - 1 - (l  + p ) , b = ( /+ 1)  Mm  , pi^  = 'o , M £ i . A plus  sign  over  a 
term  indicates  that  its  value  is  to  be  updated  on  the  current  iteration,  or  set 
equal  to  its  value  on  the  current  iteration,  if  it  has  already  been  updated, 
as  in  the  case  of  the  term  "P*.,  • The  quantity  coe  is  the  relaxation  factor 
used  at  elliptic  points. 
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Line  Relaxation 


The  sequence  of  iterations  to  be  followed  in  the  relaxation  solution 
is  as  follows:  starting  at  a given  radius,  all  the  points  on  the  line 

£=  constant  will  be  solved  for  by  line  relaxation.  This  line  will  be  swept 
downstream,  and  the  sweep  repeated  a number  of  times  before  moving  to  the  next 
radius.  Thus,  it  is  convenient  to  arrange  the  potential  equation  in  a tridiagonal 
form 


n f i-  - *»  rJ  T i-  - 1 l 

Rk  <t>K  * + C«  = Dk  . L = -2,3,...,  L MX  - / (11) 

where  LAfX  is  the  index  corresponding  to  £ = In/ 6 . The  coefficients  > 

L L.  i _ 

3*  , C*  , and  O*  are  all  known,  in  terms  of  the  independent  variables, 

and  values  of  <£  at  neighboring  points  in  the  grid  found  on  the  previous 
i-  sweep  (or,  in  the  case  of  some  terms  at  K-i  , found  on  the  current  2-  sweep). 
The  explicit  forms  of  these  coefficients  depend  on  whether  the  flow  is  locally 
hyperbolic  or  elliptic.  This  is  decided  by  the  sign  of  the  centered-difference 
form  of  the  coefficient  cf  (f>zz  : 

' * -2AT 

where  is  evaluated  from  the  values  that  were  found  on  the  previous 

2 - sweep . If  \J  is  positive,  the  point  is  elliptic,  and  centered  differences 

are  used  throughout;  if  V is  negative,  the  point  is  hyperbolic,  and  upwind 
differences  are  used  for  zl  • 


The  tridiagonal  solution  is  found  in  the  usual  way  (see  Reference  15, 
for  example) . Let 


v - 

T K 

£ (L)  Nd>‘-+'  + F (L) 

T K 

(13) 

Then  the  quantities  £ and  F obey  the  recursion  relations: 

E(L)  - 

-K/iK  + cLK£U-n] 

(14) 

F(L)  = 

IX  + CLK£(L-t)] 

(15) 
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The  boundary  or  periodicity  conditions  at  L - / are  used  to  set  E O)  and 
F <»)  . Then  the  recursion  formulas  [Equations  (14)  - (15)]  give  £ (.l)  and 
F (L)  for  L = 2 , 3 lux  - t . Boundary  or  periodicity  conditions  at  L = LM x 

are  then  used  to  set  (f>  at  L MX  , and  then  Equation  (13)  is  used  to  work  back 
down  to  L -2  . After  these  (tentative)  values  of  0 have  been  found,  (and 
stored  in  an  array  5(0  ) the  solution  is  updated  by  the  relaxation  formula 

"0 L = ujS(L)  + 0-uO  y01'  (16) 

where  the  values  used  for  the  last  term  are  those  from  the  previous  iteration. 

The  next  step  is  to  use  the  periodicity  conditions  (upstream  and 
downstream  of  the  blades)  and  the  blade-surface  boundary  conditions  (within 
the  blade  row)  in  order  to  set  the  proper  values  at  1 = 1 and  L = LHX . These 
are  taken  up,  in  turn,  in  the  paragraphs  below. 


Periodicity  Conditions 

1.  Upstream  of  the  blade  row: 


(a)  here  E(l)  is  set  equal  to  zero,  and  F (1 ) is  set  equal  to  (j> 


(b)  at  £ - 2tt/3  : here  the  proper  way  to  enforce  periodicity 


is  to  re-write  Equation  (11) , with  0 


A/  j - > 


replaced 


, v , 2 
by 


LHX  vt  ‘-HX.  « 1 i.MX-1  n4-"*  N . * fl  71 

BX  <t>K  + CK  <Pk  = DX  ~ RX  4*  U J 


V,  2 


The  quantity  # 


LH*  -1 


is  then  eliminated  in  favor  of 


'a.  *■' 

K 


by  using  Equation  (13),  resulting  in  an 


explicit  relation  for  v0 


I .MX 

k • 


2.  Downstream  of  the  blade  row: 

In  this  region,  it  is  necessary  to  account  for  the  jump  in 
potential  and  its  radial  derivative  across  the  trailing  vortex  sheet  (see 
Reference  8) . 


10 


I 


(a)  here  E(l)  is  set  equal  to  zero,  arul  F(l)  is  set  equal  to 

(b)  at  2?  - 2n/B  : in  order  to  evaluate  the  even  part  of  the 

'C,  -derivatives  correctly,  it  is  necessary  to  replace 
"W"  by  < ~ A $ (see  Figure  4).  The  odd  part 
also  makes  a contribution  to  <j> (see  Equation  (54), 
Reference  8);  the  net  result  is 


« ' jo+f1) 


A cj>(AJ  + 0 -2AiJ(W)+A<J(/V-0 
(A/O)2 


a4,(n  + i)-a4(n-,)\  ' Ul  + 

2j0  Ap  J (A%)X 

This  formula  is  then  used  to  derive  an  equation  similar  to 

V t <-MX 

Equation  (17),  from  which  an  explicit  relation  for  <t> K 
can  be  found. 


Blade-Surface  Boundary  Conditions 


There  are  two  options  for  specifying  boundary  conditions  at  the  blade 
surfaces,  designated  by  Roman  numerals  I and  II: 

I.  Here  the  blade  shape  is  specified.  Thus,  the  surface  slope 
(see  Reference  8,  Equation  (43))  is  given  at  all  grid  points  on  the  suction  and 
pressure  sides  of  the  blades. 


(a)  at  <C  = o : 

±L  \ 

dsJ~ 


ji_  \ _L  JLt.  

^ a ^ L * 1 P 3%  » *P X 


1 _34_ 

px  32 


Thus , 

JUI  clri  j P2  H 

2Z  L _ ,0  d. s L + i+p2  n „ 

£ =0  h~o  r c;  = o 

This  is  used  to  approximate  , as  follows: 

* l “k , , 1 


oLrj 
O 

' d s 


i 
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When  this  is  substituted  into  the  potential  equation,  at 
2f-o  , and  the  result  compared  with  Equation  (13),  an 
explicit  relation  for  w cfr  'K  results. 

(b)  at  % - zn  J 8 : an  analogous  procedure  is  used  here,  the 

major  difference  being  that  the  if  -derivatives  are 
replaced  by: 

, I - — + ^ 1 **  ^ j 

* G - ^ i + p2 

Thus, 

mIlh*  _ N 1 LM*  " 1 . 

IF1 — ) (20) 

The  analogous  sequence  of  steps  then  leads  to  an  expression 
for 


II.  Here  the  given  quantities  are  values,  at  all  grid  points  on  the 
blade  surfaces,  of  the  loading  ACP  and  the  thickness  t , defined  as:  (see 
Reference  8,  Equation  (43)) 


A Cp 


+L-ru 

f 0O  &OO 


(21) 


= -2  (l  + p*) 


4>: 


i 1 4-M  * 


1 + 


P‘ 


-*4>, 


+ z<t>. 


L = t 


t'  s>  = in-  _ 

d 5 d 5 

The  first  of  these  is  used  at  L - 


‘An 

W . 


w „ 


, the  second  atLMx: 


(22) 


(a) 


at  if  -o  : Equation  (21)  can  be  rewritten  as: 


rJ  , I W,  I 

<PK  - <PLE  ~ V*  ~V,_B 
But,  from  periodicity  conditions 
the  potential  at  ? = o can  be  set  from  the  formula: 
‘‘2  ACP 


bf . tnx  N.  LMX 

- <*.«  + 


oL  Z 


A(2P 

2 

"P  ’ 

r LI? 


AU  i-hx 
<PlE 


so 


w 1 1 


9, 


* , LHX 

Q>  +• 


d l 


(23) 


The  integral  appearing  here  is  done  once,  at  the  start  of 
the  calculations. 
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(b)  at  %'2n/8  : Equation  (22)  is  rearranged,  by  expressing 

u„  in  terms  of  and  ; then  <pi  is  replaced  in 

terms  of  A CP  . The  result  is 

= <kL.,  - ' ft  C5)  c24) 

This  expression  is  then  used  to  approximate  . as 

in  Equation  (20),  and  the  result  is  substituted  into  the 
potential  equation  at  LMx,  leading  to  an  explicit  relation 

C V T LH* 

for  cj>  < 

Blade  Loading 

The  force  per  unit  span  acting  normal  to  the  blade  chord  can  be 
related  to  the  jump  in  potential  (see  Figure  5). 

~~  * f (j>u  - fu.)  dS  = 2 T foo  aoo  V'  + P*  ' 77-  A<p(/°>  (25) 


The  quantity  A<p  is  also  related  to  the  circulation 
r _ f f..  I ..1  \ 1 c a<*> 


r = f KL  - a*L) d5 


Blade  Shape 

When  the  calculations  are  done  with  prescribed  thickness  and  loading 
distributions,  the  blade  shape  must  be  found  after  the  solution  converges. 
This  is  done  by  integrating  the  condition 

oM-  “-n  Oul.l  lx- 00  r r f \ . , 


V — ) 

J o ljl,  L 


The  trapezoidal  rule  was  used  for  these  integrations. 


Far-Field  Conditions 

The  boundary  condition  used  at  £—►-«»  is  that  the  perturbation  po- 
tential be  zero.  This  is  done  by  setting  "(p^Oior  all  N and  L , where 
the  station  corresponding  to  K * f is  taken  several  chord  lengths  upstream 
of  the  leading  edge. 
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The  condition  that  there  be  no  disturbances  upstream  is  correct  for 
flows  where  u/0 /a^  < 1 • However,  for  the  supersonic-tip  case,  where  VJ0  /a.„  > l, 
the  proper  condition  is  a radiation  condition,  i.e.,  for  c » , it  is 

required  that  any  waves  present  in  the  flow  must  have  come  from  the  blade  row. 
At  the  present  time,  this  condition  has  not  been  incorporated;  thus,  the 
solutions  are  restricted  to  cases  where  tip  Mach  number  is  subsonic. 


At  large  distances  downstream  of  the  blade  row,  the  perturbation 
quantities  are  required  to  be  independent  of  £ , at  constant  £ , i.e.,  they 

do  not  change  in  the  direction  along  the  helical  paths.  This  is  enforced  by 
selecting  values  of  cp  at  K - KHX,  equal  to  the  constant  value  required  for 
mass  conservation  (see  Reference  8,  Equation  (56)). 


where 


-k. 


DRF 


V 

T K I 


-DBF-  v4>; 


KMX  - 1 


KMX 


DCF 


(28) 


DRF 

DBF 

DCF 


KMX  ^ KMX  - KMX  ~ ^ KMX  -X  ) 


(2-  KMX  ~ KMX  -i)  / (**  KMX  - 1 '^KMX  - KMX  -I  ^~KMK  ) 


( KMX  £ KMX  -1  ~ ^ KMX  -X  )/  ^ KMX  ~ ^ KM  X - 1 ) 6^"  KMX  ~ £ KMX  -1  ) 


Hub  and  Shroud  Conditions 

In  general,  the  boundary  conditions  at  the  hub  and  shroud  require 
that  the  radial  derivative  of  the  potential  be  proportional  to  the  slope  of 
the  surface.  For  the  present  formulation,  these  surfaces  are  taken  to  be 
cylinders;  thus,  the  requirement  is  that  d<P/dp  = 0 . This  is  enforced  by 
setting  <p  at  N-t  and  A/*/vmx  such  that  the  derivative  is  zero,  i.e., 


4 

3 


V, 


KMX  L 

<P 

^ K 


4 KMX-1  L 

1 V* 


f KMX-1  i- 

1 


(29) 
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Updating  of  the  Jump  in  Potential 


When  the  calculations  are  being  done  in  the  mode  where  the  geometry 
is  given,  the  quantity  Acf>(/b)  is  updated  as  the  calculation  proceeds,  in  the 

,(16) 

+ tJCa,  in 

- > p ’ sr 


manner  outlined  by  Ballhaus  and  Bailey 

UJC*.  \ T / 

V ' *( 


A 4>  (f)  = uJf 


)v  Dai  i 

{*( 


u. 


P > 


t ( i - u>p ) a 4>  cp) 


The  relaxation  factor  u)P  used  here  is  taken  equal  to  3/kmx,  since  a monotone 
convergent  sequence  of  iterations  will  behave  according  to 


- uJ P • I TK.- 


V 


where  ITK  counts  the  iterations  in  the  axial  direction.  Since  the  total 
number  of  sweeps  in  the  axial  directions  is  usually  chosen  to  be  of  the  order 
of  KMx.  , the  above  choice  of  uJP  will  assure  that  the  argument  of  the  expo- 
nential factor  is  sufficiently  large. 


(30) 


Because  the  perturbation  velocity  components  are  forced  to  be  periodic 
downstream  of  the  trailing  edge,  the  values  of  <P  at  the  trailing  edge  that 
enter  Eq.  30  will  tend  toward  the  value  of  A 4>  required  by  the  Kutta  condition, 
i.e.,  the  value  that  makes  the  perturbation  pressures  at  and  £ = m//5  equal 

This  condition  can  be  enforced  more  strongly  by  selecting  A , at  each 
iteration,  to  be  precisely  the  value  that  will  make  the  perturbation  pressures 
equal  at  the  trailing  edge.  If  subscripts  a and  b denote  the  locations  of  the 
last  two  points  on  the  airfoil,  this  value  is: 

x*  i~(2TE  ^4>a.  +(2re 
A Q>  = A 

( ~ ) ( 2 2-TE  ~ ol~  2- b) 

where  ZT£  - u>C A/  a m 

= & O-x.b’P'0')  ~ <t>(*cL,b  » P>  Jf-) 

The  computer  program  contains  an  option  for  using  this  formula.  Under  this 
option,  60p  is  read  in,  and  A <p  is  updated  by 

*■ 

A<f>  = ^ <*>  * ( ' - tJp)  ^ *+> 
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A limited  number  of  runs,  made  with  this  option,  showed  only  minor  differences 
from  those  which  used  Eq.  30.  All  of  the  results  presented  below  used  Eq.  30. 


Grid  Used 

The  grid  used  for  the  calculations  reported  here  had  a uniform  spacing 
in  the  p-  and  t,  -directions.  The  grid  was  unevenly  spaced  in  the  2 -direction, 
as  follows:  the  upstream  and  downstream  boundaries  were  set  at 


-2  

^ CO 


Then  the  intermediate  points  were  located  at 


where 


i(K)  = ZM  + 


— JU, 
20c 


K A Z 
2 - KA  X 


< = 2}..-,  KMX. 


J (*z 


At  - 2/  (K  MX  + t ) 


2 Oc  = 


JU  [A  Z / (2  - At  )] 


Relaxation  Factors 

In  updating  the  points  on  a given  line  ( 2 = constant,  p = constant, 
O * £ 5 iv./ B ) the  relaxation  factor  was  taken  as  either  of  two  values, 
depending  on  whether  there  were  any  hyperbolic  points  on  that  line.  If  there 
were,  a hyperbolic  relaxation  factor  was  used;  if  all  the  points  were  elliptic, 
an  elliptic  relaxation  factor  was  used. 


I 
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The  hyperbolic  relaxation  factor,  called  RxH  , was  read  in  as  a 
constant,  less  than  1.0.  The  elliptic  relaxation  factor  was  taken  as  a 
function  of  2 , which  varied  from  an  input  maximum  value,  called  RxE  , t.o 

the  value  1 . 0 at  large  distances  upstream  and  downstream  of  the  blades, 
according  to  the  formula: 

r /z-  -f  \4 ) 

- i-°  + (RXE  - t)  i h (32) 

[ \ uyC^ /tf„  / J 

When  this  was  not  done,  it  was  found  that  values  of  the  potential  far  upstream 
and  far  downstream  did  not  converge. 


II.  RESULTS 


The  computer  program  described  above  has  been  used  to  carry  out  two 
calculations,  for  the  purpose  of  demonstrating  its  capabilities.*  The  blade  row 
used  for  both  of  these  calculations  is  shown  in  Figure  6.  It  had  30  blades, 
with  a hub-to-tip  radius  ratio  of  0.5  and  a solidity  Co,/  LT  of  0.5.  Both 
calculations  used  an  axial  Mach  number  of  0.4,  and  an  angular  Mach  number  at 
the  tip  of  0.8,  making  the  resultant  Mach  number  at  the  tip  0.894.  The  grid 
used  for  these  calculations  had  40  points  in  the  1 -direction,  20  points  in 
the  P -direction,  and  10  points  radially,  for  a total  of  8000  grid  points. 
Constant  grid  spacings  were  used  in  the  p - and  £,  -directions;  a variable  grid 
was  used  in  the  2 -direction,  with  points  more  closely  spaced  in  the  vicinity 
of  the  blades  (14  of  the  40  points  lay  between  the  leading-  and  trailing-edge 
stations).  The  grid  extended  from  2 Co.  upstream  of  the  leading-edge  station 
to ZCa.  downstream  of  the  trailing-edge  station,  as  noted  in  Section  I. 

The  relaxation  procedure  began  by  finding  the  solution  at  the  hub. 

Line  relaxation  was  used,  with  all  points  on  the  line  2 = constant,  p = con- 
stant updated  simultaneously.  This  line  is  then  swept  from  upstream  to  down- 
stream, and  the  sweep  is  repeated  a number  of  times,  typically  the  same  as  the 
number  of  grid  points  in  the  £ -direction,  in  order  for  information  to  be 
carried  from  one  end  of  the  grid  to  the  other.  The  solution  then  proceeds  to 
the  next  radius,  where  the  process  is  repeated.  After  this  first  sequence  of 
radial  iterations  has  reached  the  tip  radius,  it  is  then  repeated.  The  number 
of  repetitions  is  typically  the  same  as  the  number  of  radial  surfaces  in  the 
grid.  This  procedure  is  analogous  to  one  which  has  been  used  for  isolated- 
wing  calculations . 17 

The  program  was  written  in  FORTRAN  IV  and,  using  the  H compiler,  was 
run  on  the  IBM  370/168  of  the  Martin-Marietta  Corporation  Great  Lakes  Data 
Center.  It  was  found  that  the  computing  time  required  was  36  microseconds  per 

k 

An  earlier  version  of  these  results,  calculated  with  a non-conservative 
difference  scheme,  was  presented  in  Ref.  16. 
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grid  point  per  visitation,  exclusive  of  the  time  required  for  output.  Each  of 
the  two  calculations  reported  below  took  about  four  minutes,  plus  the  time  used 
in  output.  (The  latter  time  requires  an  additional  5 to  10  seconds,  depending 
on  the  amount  of  output  desired.) 


First  Case:  Given  Blade  Shape 

In  the  first  demonstration  calculation,  the  blade  geometry  was 
specified.  The  blades  were  taken  to  have  parabolic-arc  distributions  of  thick- 
ness and  camber  (see  Figure  7). 


lu.  (5) 

L 

• * cs>  1 i 

t < s ) 

t (S)  = 

4 ^ majc  if'’)  | 

| S t^'5]} 

[c  <r> J * 1 

-A(S)  = 

[c<rj]3 

[s  [c-s] j 

The  maximum  thickness  was  chosen  as  a constant,  of  such  a magnitude  that  the 
variation  of  chord  length  with  radius  produces  a thickness-to-chord  ratio 
tnxxx/Cc.r'>  which  varies  from  6%  at  the  tip  to  9.49%  at  the  hub.  The  maximum 
camber  was  also  chosen  as  a constant,  of  such  a magnitude  as  to  make  the  camber 
vary  from  4%  of  the  chord  at  the  tip  to  6.33%  of  the  chord  at  the  hub.  The 
blade  shapes  at  the  hub  and  tip  are  shown  in  Figure  8. 

The  convergence  of  the  solution  is  shown  in  Figure  9,  which  presents 
the  variation  of  the  streamwise  perturbation  velocity  at  mid-chord,  on  the 
suction  side  of  the  blade.  The  iteration  number  shown  here  counts  the  number 
of  times  the  values  at  a given  point  have  been  updated.  The  calculation  was 
done  in  ten  successive  stages:  most  of  these  used  two  radial  sweeps,  with 
20  axial  sweeps  on  each  radial  one;  on  several  of  the  stages,  only  one  radial 
sweep  was  used,  with  40  axial  sweeps. 

The  relaxation  factor  (see  Equation  (10)  was  taken  as  1.0;  the 
quantity  Rxh  was  chosen  as  0.9,  while  ff-XE  was  set  equal  to  1.0  for  the  first 
120  iterations,  and  equal  to  1.2  thereafter.  Several  attempts  were  made,  at 
later  stages  of  the  solution,  to  increase  RxL  to  1.4,  but  these  led  to  divergent 
behavior. 
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The  results  of  this  calculation  show  that  the  flow  accelerates  to 
locally  supersonic  conditions  over  the  outer  half  of  the  span.  Contours  of 
the  local  Mach  number  are  shown  in  Figures  10a  and  10b  for  radial  stations  near 
the  hub  and  near  the  tip. 


Second  Case:  Given  Thickness  and  Loading 

The  second  demonstration  calculation  was  done  using  the  computer  pro- 
gram in  the  design  mode.  The  thickness  distribution  was  chosen  to  be  the  same 
as  used  in  the  above  calculation,  while  the  loading  distribution  had  the  form 
shown  in  Figure  11.  The  value  of  A Cpo  if)  varied  linearly  between  pre- 
scribed values  at  the  hub  and  tip.*  The  values  used  were  &.  = 0.6,  ACpa  = 0.5 
at  the  hub,  and  1.0  at  the  tip.  The  flow  resulting  from  these  specifications 
was  everywhere  subsonic;  contours  of  the  local  Mach  number  at  two  radial 
stations  are  given  in  Figures  12a  and  12b. 


The  convergence  of  the  solution  is  shown  in  Figure  13,  using  the  same 
location  as  in  Figure  9.  The  relaxation  factors  u)e  and  RxH  were  again  chosen 
as  1.0  and  0.9.  However,  it  was  found  that  the  relaxation  factor  used  for 
elliptic  points  had  to  be  quite  small  at  the  beginning  of  the  calculation.  The 
values  used  were: 


Iteration  Number 

RXb 

1-200 

0.1 

201-400 

0.4 

401-1000 

0.8 

The  blade  shape  required  to  achieve  this  loading  is  shown  in 
Figure  14  for  the  hub  and  tip.  It  is  interesting  to  note  that  an  angle  of 
attack  makes  up  the  largest  portion  of  the  blade  shape,  with  only  a modest 
contribution  from  camber. 


Strictly  speaking,  this  distribution  is  inconsistent  with  the  condition  v * o 
at  the  hub  and  tip,  since  it  is  always  true  that 
v( Z , p , o)  - v (i , p,  Ad>/3 

As  a consequence,  the  radial  velocities  Calculated  for  this  case  are  inaccurate 
near  the  hub  and  tip. 
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III.  CONCLUDING  REMARKS 


The  principal  significance  of  the  results  presented  above  is  that  the 
first  step  spoken  of  in  the  introduction  has  been  carried  out:  it  is  now 

possible  to  make  three-dimensional  transonic  flow  calculations  of  periodic 
flows  in  a compressor  blade  row.  Because  a relaxation  method  is  used,  the 
computing  times  required  are  quite  modest.  Further  steps  are  now  possible,  in 
the  directions  of  application  and  extension. 

Some  of  the  applications  that  need  to  be  made  are  to  compare  the  pre- 
dictions of  this  method  with  experiment.  As  noted  earlier,  it  can  be  expected 
that  the  sma 1 1 -disturbance  theory  will  have  a certain  range  of  validity;  the 
extent  of  this  range,  and  the  principal  factors  which  limit  it  remain  to  be 
determined  by  comparisons  with  experimental  data.  In  addition,  there  are 
several  questions  that  can  be  answered  by  further  application:  among  these 

are  the  effect  of  blade  shape  and  loading  on  the  location  of  supersonic  zones, 
the  types  of  camber- line  distributions  required  to  achieve  various  pressure 
distributions,  and  the  degree  of  validity  of  quasi-two-dimensional  solution 
methods.  All  of  these  questions  can  now  be  answered,  in  the  small-disturbance 
approximation,  by  conducting  a series  of  calculations  in  which  the  appropriate 
parameters  are  varied.  The  experience  gained  to  date  suggests  that  these 
calculations  can  be  done  at  quite  reasonable  cost. 

In  addition,  the  availability  of  a three-dimensional  inviscid  solution 
makes  it  possible  to  consider  the  three-dimensional  boundary- layer  flow  that 
would  result.  By  examining  the  pressure  distributions  on  the  hub,  shroud,  and 
blade  surfaces,  the  general  features  of  these  boundary  layers  can  be  described, 
and  quantitative  estimates  of  viscous  losses  can  be  made. 

Looking  beyond  these  applications , there  are  many  extensions  that  can 
be  made.  Some  of  the  most  important  ones  are  the  inclusion  of  further  non- 
linear terms  in  the  potential  equation,  which  may  enable  consideration  of 
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highly  cambered  blades,  large  turning  angles,  and  pressure  ratios  substantial ly 
larger  than  1.0.  These  extensions  might  also  allow  incorporation  of  further 
geometrical  complications,  such  as  rounded  leading  edges,  part-span  shrouds, 
tapered  blades,  and  variations  in  the  annular  cross-section. 

One  of  the  major  limitations  of  the  present  results  is  its  poor 

representation  of  shock  transitions.  This  problem  can  be  overcome  by  the  addi- 

18  19 

tion  of  a shock-fitting  technique.  ’ Once  that  has  been  done,  the  necessary 
tools  will  be  at  hand  for  studies  of  shock-wave/boundary  layer  interaction  in  a 
transonic  compressor  rotor. 

Even  with  the  full  nonlinear  potential  equation,  there  are  still  a 
number  of  limiting  assumptions,  especially  in  regard  to  the  absence  of  gradients 
in  the  inlet  flow.  An  alternate  formulation,  recently  been  developed  by 
McCune  and  Hawthorne,  is  free  from  this  restriction.  This  formulation 

uses  small  perturbationsof  the  axisymmetric  mean  flow;  it  may  be  that  the 
numerical  solution  of  the  perturbation  equations  in  this  formulation  would  be 
a profitable  avenue  to  follow. 

On  an  even  longer  time  scale,  there  are  many  other  techniques  that 
have  proven  to  be  of  great  value  in  relaxation  calculations  of  external  flows, 
and  which  hold  comparable  promise  for  internal  flows.  For  example,  the  area 
of  convergence-acceleration  techniques can  now  be  considered  for  three- 

2 S 27 

dimensional  turbomachinery  flows,  and  use  of  a velocity-component  formulation, 
capable  of  treating  flows  with  radial  inlet  gradients,  can  be  contemplated. 

As  was  the  case  with  external-flow  calculations,  the  experience  gained  in  the 
small-disturbance  approximation,  in  regard  to  such  matters  as  coordinate 
systems,  periodicity,  and  rotated  differences,  helps  to  show  how  further  ad- 
vances can  be  made. 
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Figure  1 BLADE-FIXED  COORDINATES 


ROTOR  IS  STATIONARY  IN  A 
HELICAL  APPROACH  FLOW 


Figure  2 COORDINATE  SYSTEM  AND  FINITE-DIFFERENCE  GRID,  SHOWING  FORWARD 
MACH  CONE  AT  POINT  K,  L.  SURFACE  SHOWN  IS  AT  CONSTANT  RADIUS 
(N  = CONSTANT) 


Figure  3 RECTANGULAR  GRID,  SHOWING  MISALIGNMENT  BETWEEN  FLOW  DIRECTION 
AND  COORDINATE  DIRECTIONS. 
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Figure  6 BLADE  ROW  USED  FOR  DEMONSTRATION  CALCULATIONS 
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Figure  7 


DEFINITIONS  OF  BLADE-SURFACE  GEOMETRY 
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